Understanding the impact of modiolus porosity on stimulation of spiral ganglion neurons by cochlear implants

Moderate-to-profound sensorineural hearing loss in humans is treatable by electrically stimulating the auditory nerve (AN) with a cochlear implant (CI). In the cochlea, the modiolus presents a porous bony interface between the CI electrode and the AN. New bone growth caused by the presence of the CI electrode or neural degeneration inflicted by ageing or otological diseases might change the effective porosity of the modiolus and, thereby, alter its electrical material properties. Using a volume conductor description of the cochlea, with the aid of a ‘mapped conductivity’ method and an ad-hoc ‘regionally kinetic’ equation system, we show that even a slight variation in modiolus porosity or pore distribution can disproportionately affect AN stimulation. Hence, because of porosity changes, an inconsistent CI performance might occur if neural degeneration or new bone growth progress after implantation. Appropriate electrical material properties in accordance with modiolar morphology and pathology should be considered in patient-specific studies. The present first-of-its-kind in-silico study advocates for contextual experimental studies to further explore the utility of modiolus porous morphology in optimising the CI outcome.


Supplementary material
A. Mesh selection for the 3D computational domain of the cochlea.
In the present study, the 3D spiral computational domain of the cochlea, with a few millimeters in length, contains thin membrane structures and as small as 2-µm axonal initial segments attached to the spiral ganglion neuron (SGN) cell body of 30-µm diameter.For the finite-element discretization, these three orders of magnitude in variations of the model geometry demand a large number of mesh elements to restore the shape of the extremely small subdomains.COMSOL Multiphysics ® (briefly called COMSOL hereafter) finite-element software offers automatic inbuilt meshing options suitable for the computational domain according to the chosen study module; in the present case, the AC/DC module.
We studied all seven automatic meshing options and related mesh parameters viz.maximum element size (Element size) that limits how large each element could be, maximum element growth rate (Growth rate) that limits the size difference of two adjacent elements, curvature factor (Curvature factor) that limits how large a mesh element can be along a curved boundary, and resolution of narrow regions (Resolution of regions) that controls the number of layers of mesh elements in a narrow region chosen by COMSOL.Supplementary Table 1 shows the above-mentioned mesh parameters, solution time, and degrees of freedom (DoF) solved for each mesh option.Taking the numerical solution (transmembrane potential of an SGN) obtained by the finest mesh (mesh number 7) as the best solution (∅  ), the relative (  ) and the absolute (  ) error magnitudes were calculated, respectively, according to equations (1) and (2), using numerical solutions obtained with the remaining meshes (∅  ).As an eighth mesh, a manual meshing option was chosen to reduce the number of mesh elements while taking care of the thinnest subdomains, the computational cost, and the accuracy of the solution.
The dependence of   and   on the mesh type is shown in Supplementary Fig. 1a.
Since the heterogeneous electric conductivity distribution was modeled on the modiolus subdomain, an adaptive mesh refinement was implemented on the modiolus to minimize the interpolation of conductivity values.Supplementary Figure 1b shows the final tetrahedral mesh obtained by the manual meshing (mesh number 8) used for all simulations in the present study.Supplementary Figure 1 | Finite-element mesh selection for the 3D computational domain.a, Relative error percentage and absolute error magnitude in numerical solutions while using automatic meshes generated by COMSOL (mesh numbers 1-7) and a manual mesh (mesh number 8).b, Finite-element mesh (mesh number 8) generated by manual meshing on the 3D computational domain.Several thin subdomains and their respective meshes are shown in the zoomed pictures.

B. Parametric study of 'regionally kinetic' porosity equations on a two-dimensional domain.
The present study proposes equations (3) and (4) to model the random pore distribution on the modiolus.This system of reaction-diffusion equations forms regionally kinetic (RK) random patterns for each time step after the spiral breakup while keeping the area occupied by the higher values of the state variable u almost constant in the computational domain.
In the present study, we assigned a = 0.64, b = 0.02, c = 0.08, D = 1, and β = 0.1.Here, we demonstrate how the newly introduced term −βu and the value k = 9 affect the system dynamics compared with the Bär and Eiswirth (BE) model and how the parameter β impacts the dynamics of the RK porosity equations on a square domain of length L = 3 mm and equidistant grid size h = 50 µm.The time step is similar to that in Barkley's model 2 ; i.e., proportional to the square of the grid size such that ∆ = 0.4ℎ 2 = 1 ns.Studying the BE and RK equation systems by choosing such a time step facilitates a better understanding of critical transition points, such as spiral wave formation, wave breakup, and random pattern formation.
The number of solved DoF was 30,000 using a time-dependent solver; the solution time for both BE and RK equations was 2.5 minutes to complete 500-time steps.
In the present study, the instants of interest are those generating random patterns by the RK equations from about the 100 th time step (0.1 µs) onwards.Supplementary Figure 2d shows the formation of random patterns for β = 0.1, which are regionally kinetic and result in almost the same area occupied by the higher values of u.Supplementary Figure 3a shows the percentage of area occupied for each indicated value of u at each time step.Supplementary Figure 3b shows the impact of β on the spread of higher values of u and subsequent increment in its occupying area in the computational domain at the 100 th time step.The β values higher than 0.2 resulted in an un-wanted spread of higher values of u (Supplementary Figure 3c).

For β = 0
, k = 6.75, and D = 1, equations (3) and (4) become the system of reactiondiffusion equations proposed by Bär & Eiswirth 1 , who have conducted several numerical experiments by choosing different values for the system parameters a, b, and c.They modified the piecewise linearized FitzHugh-Nagumo equations proposed byBarkley 2 to study spiralbreakup and turbulence patterns in an excitable media.Bär & Eiswirth studied the reactiondiffusion equation system proposed in 1 in a two-dimensional computational domain to reduce the computational cost while providing essential insights into the system dynamics.The parameter values tested through such studies can be used for the three-dimensional domain.
Hence, β should be chosen between 0.1 and 0.2, inclusive, for the RK equation system in the present application.Based on insights acquired from this comparative study, we assigned a = 0.64, b = 0.02, c = 0.08, D = 1, and β = 0.1 to solve the RK equations with the time step of 0.1 µs on the 3D computational domain of the modiolus.Supplementary Figure 3 | System properties of RK equations.a, Almost constant percentage of area (% u) occupied by the indicated higher values of the state parameter u at each time step (t) given in Supplementary Figure 2d.b & c, For t = 0.1 µs (100 th time step), with small increments of β between 0.1 and 0.2, the higher values of u occupied more area in the computational domain.β = 0.3 resulted in an unwanted spread of u values.